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ABSTRACT 

We  investigate  the  possibility  of  exploiting  the  properties  of  a  detected  Low 
Probability  of  Intercept  (LPI)  signal  waveform  to  estimate  time  delay,  and  by 
geometry,  angle  of  arrival.  We  consider  the  case  where  a  highly  correlated 
signal  is  received  at  two  stationary  passive  receivers.  The  signal  source  is 
assumed  stationary,  and  the  signal  waveform  designed  so  that  the  ambiguity 
function  has  a  sharp  peak.  We  also  restrict  ourselves  to  low  signal-to-noise 
ratios,  namely  10  dB  and  less.  We  also  examine  the  minimum  time-delay 
estimate  error  -  the  Cramer-Rao  bound. 

The  results  indicate  that  the  method  works  well  for  highly  correlated  pulsed 
signals,  and  may  prove  useful  for  other  types  of  signals,  such  as  CW  signals 
and  pseudo-random  noise. 
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Time  Delay  Estimation 


EXECUTIVE  SUMMARY 

The  theory  of  radar  signal  (target  echo)  detection  can  be  formulated  with  the  use  of 
the  ambiguity  function.  Such  functions  are  designed  to  have  sharp  peaks,  and  a  high 
correlation  on  return.  Resolution  is  determined  by  the  width  of  the  peak. 

To  improve  detection  and  target  resolution  such  signals  make  use  of  pulse  compression 
techniques. 

An  Electronic  Support  measures  (ESM)  system  can  exploit  this  information  by  using  two 
passive  receivers,  separated  in  space.  Knowing  the  received  waveform  will  have  a  sharply 
peaked  ambiguity  function,  we  look  at  the  correlation  between  the  output  of  one  receiver, 
with  the  output  of  the  second. 

This  correlation  will  give  us  the  time-delay,  to  the  same  resolution  as  was  “built  into”  the 
original  waveform.  From  this,  using  simple  geometry,  the  signals  angle  of  arrival  can  be 
determined. 

Waveform  design  is  performed  in  the  frequency  domain,  making  use  of  complex  signals. 
As  a  result,  we  work  primarily  within  the  frequency  domain. 

The  starting  point  is  the  analytic  signal.  We  generate  such  a  signal,  and  form  the  product 
of  the  signal  with  its  (complex  conjugate)  time-delayed  version.  This  product,  in  the 
frequency  domain,  is  the  cross  spectral  density.  Transforming  to  the  time  domain,  we  find 
the  peak  of  the  time-domain  correlation  function,  using  quadratic  interpolation. 

This  peak  gives  the  estimate  of  the  time-delay. 

We  have  run  simulations  over  various  signal-to-noise  ratios  (SNRs)  with  continuous  wave 
and  single  pulse  chirped  signals.  As  we  are  interested  in  the  ability  to  estimate  time-delay 
at  low  SNR,  we  run  our  simulations  with  SNR  at  10  dB  and  less. 

The  continuous  wave  case  allows  us  to  observe  the  effects  of  frequency  wrap  around,  as 
the  waveform  will  resemble  a  partial  sawtooth  during  the  sample  time. 

The  continuous  wave  simulations  show  that  there  is  a  bias  due  to  the  frequency  wrap 
around,  producing  a  rough  measure  of  the  actual  time-delay.  The  single  pulse  case  gives 
very  good  estimates  of  the  time  delay,  although  all  estimates  degrade  as  the  SNR  decreases. 

We  have  also  tested  the  algorithm  on  pseudo-random  noise  generated  during  the  Electronic 
Warfare  and  Radar  Division  (EWRD)  Durex  trial.  Again,  although  we  do  not  accurately 
estimate  time-delay,  we  do  get  a  rough  measure  of  the  delay. 

The  results  indicate  that  the  method  works  well  for  highly  correlated  pulsed  signals,  and 
may  prove  useful  for  other  types  of  signals,  such  as  CW  signals  and  pseudo-random  noise. 
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1  Introduction 

The  ability  to  find  the  direction  of  an  emitted  signal  is  of  importance  in  a  wide  range  of 
fields,  including  electronic  intelligence  (ELINT)  and  electronic  warfare  support  (ES). 

For  a  digital  receiver,  direction  finding  is  necessarily  preceeded  by  signal  detection,  and 
parameter  estimation,  and  this  leads  to  the  theory  of  waveform  design.  In  developing  the 
theory  of  waveform  design  and  analysis,  it  is  found  that  a  simplification  of  the  analysis  can 
be  achieved  by  working  in  the  complex  domain,  and  by  generalising  the  phasor  description 
of  signals  [2].  As  a  result,  it  is  sensible  to  work  in  the  frequency  domain,  returning  to  the 
time  domain  as  needed. 

For  ranging,  an  emitted  waveform  is  correlated  with  the  return  from  a  target.  As  the 
waveform  used  is  designed  to  be  highly  correlated  -  its  ambiguity  function  will  have  a 
sharp  peak  -  there  will  be  a  sharp  peak  in  the  correlation  function.  This  peak  gives  the 
time-delay,  or  range  of  the  target. 

An  electronic  support  measures  (ESM)  system  can  exploit  this  information,  and  use  it  to 
determine  the  direction  of  the  signal  source.  We  do  this  by  correlating  the  received  signal 
in  a  passive  receiver,  with  the  time-delayed  signal  in  a  second  passive  receiver,  located  a 
short  distance  away  in  space.  The  output  of  this  correlator  receiver  is  the  time-delay,  and 
by  geometry,  the  angle  of  arrival. 

This  is  shown  in  Chapter  2.  We  start  by  introducing  the  ambiguity  function,  and  the 
method  used  for  estimating  the  time-delay.  We  run  simulations  over  various  Signal-to- 
Noise  ratios  (SNR),  for  both  continuous  wave  (CW),  and  single-pulse  signals.  In  this 
chapter  we  also  derive  expressions  for  the  Cramer-Rao  bound  of  the  variance  of  the  esti¬ 
mates. 

In  Chapter  3  we  apply  the  techniques  developed  in  Chapter  2  to  pseudo-random  noise, 
generated  at  the  Electronic  Warfare  and  Radar  Division  (EWRD)  Durex  trial  [1]. 

In  Chapter  4,  we  briefly  review  a  technique  which  uses  the  phase  slope  of  the  cross  spectral 
density  to  estimate  time  delay.  In  this  way,  we  see  that  it  is  not  necessary  to  transform 
back  to  the  time-domain  in  order  to  estimate  the  time-delay. 

In  Appendix  A  we  review  an  alternative  method  for  deriving  the  Cramer-Rao  bound  for 
the  estimate  of  angle  of  arrival,  and  review  the  generalisation  of  this  method  in  Appendix 
B. 
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2  Time— Delay  Estimation 

2.1  Introduction 

Suppose  our  aim  is  to  distinguish  different  ranges  (time-delays)  at  a  receiver.  The  trans¬ 
mitted  signal  waveform  must  have  the  property  that  it  is  as  different  from  its  shifted  self 
as  possible  [3]. 

That  is,  if  ip(t)  represents  our  complex  transmitted  waveform,  then  the  mean  square 
difference  [4] 

J  \^(t)  ~  i/>(t  +  r)|2  dt  (1) 

must  be  as  large  as  possible  over  the  time  range  containing  the  delay  r. 

This  excludes  a  small  range  of  r  near  zero,  where  must  resemble  tjj(t  +  r). 

Expanding  the  term  inside  the  integral,  for  complex  waveforms,  we  find  that  we  need  to 
minimise  (except  near  t  =  0) 

9i  J  ip(t)  ip*  (t  +  t)  dt  (2) 

with  3?  denoting  the  real  part  of  the  integral,  and  *  denoting  complex  conjugation. 

We  see  this  by  considering  ^(t)  and  ^{t  +  r)  as  complex  vectors,  with  the  same  origin. 
We  keep  ip(t)  fixed,  and  rotate  ip(t  +  t)  until  the  vector  connecting  their  end  points, 
\ip(t)  —  ijj(t  +  r)|,  is  maximised. 

If  we  write 

^{t)  =  u(t )  eiujt  (3) 

then  we  minimise 

3?  e~luJT  J  u(t)  u*  [t  +  t)  dt  (4) 

an  oscillatory  function  of  r.  That  is,  if  this  becomes  negative  for  a  particular  value  of  r, 
then  there  will  be  a  corresponding  positive  value  close  to  it.  Thus  we  make  the  modulus 
as  small  as  possible. 

That  is,  we  minimise  the  modulus  of 


c(t)  =  J  u(t)  u*  (t  +  r)  dt  (5) 

which  is  the  complex  auto-correlation  function. 

We  choose  u{t)  such  that  |c(r)|  ~  0  except  at  r  =  0,  where  it  is  a  maximum.  Note  that, 
if  |c(t)|  =  0  except  at  r  =  0,  we  avoid  range  ambiguities. 

Range  accuracy  for  a  stationary  target  depends  only  upon  the  energy  of  the  signal  and 
the  bandwidth  of  the  signal  energy  spectrum,  denoted  by  f3  [4], 

For  a  single  target,  the  accuracy  with  which  r  can  be  determined  is  limited  only  by  the 
amount  of  noise  present.  If  u(t)  and  u(t  +  r)  differ,  the  difference  can  be  observed  if  noise 
is  small  enough. 
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Accuracy  of  measurement  depends  only  on  SNR  (assumed  large)  and  the  properties  of 
c(r)  near  r  =  0.  A  Taylor  expansion  gives  the  parabola 

c(r)  «  a  —  br2  (6) 


where  a  and  b  are  constants,  and  accuracy  is  found  to  be  [4] 


(7) 


where,  for  a  pulsed  radar,  1  / /?  oc  pulse  length. 

The  above  is  for  a  stationary  target. 

We  now  consider  a  moving  target.  Theoretically,  there  is  no  limitation  on  the  accuracy  of 
range  or  velocity  (frequency),  given  a  broad  spectrum  (for  range)  and  a  long  time  duration 
waveform  (for  velocity). 

However,  there  are  accuracy  limitations  on  the  combined  resolution  in  time  and  frequency. 

For  targets  at  different  ranges  and  velocities,  neither  of  which  is  known,  we  require  a 
correlation  function  for  a  combined  time  and  frequency  shift.  This  is 


x(t,  (/))  =  J  u(t)u*(t  +  r)e  l2n^>t  dt 


(8) 


with  normalisation 

x(0,0)  =  l  (9) 


The  interpretation  is  that  a  target  with  range  and  velocity  (to,  4>o)  cannot  be  distinguished 
from  a  target  shifted  in  time  and  frequency  at  (to  +  r,  +  <f>)  if  x(t,  4>)  =  1  (complete 
ambiguity) . 

That  is,  if  the  ambiguity  function  has  a  sharp  peak  at  the  origin,  and  no  other  ‘large’ 
peaks,  then  the  waveform  has  good  simultaneous  range  and  velocity  resolution. 

Resolving  a  pair  of  targets  is  not  as  straight  forward.  We  can  no  longer  expand  c(r),  or 
more  generally  y(r,  (j>),  as  a  parabola,  as  the  sum  of  two  parabolic  functions  is  another 
parabolic  function,  with  a  single  peak.  Thus  it  will  look  like  we  have  a  single  target. 

This  means  the  first  two  terms  of  the  Taylor  expansion  say  nothing  about  the  resolving 
power  of  the  waveform. 

Including  higher  order  terms,  we  have  [5] 

1x0 r,  <f>)\  =  x(0, 0)[1  -  ^(/?2r2  +  2 \A12T(j)  +  a2(j) 2  +  •••)]  (10) 


with  (3  the  effective  bandwidth  of  the  signal  when  the  mean  frequency  is  zero,  a  is  the 
effective  duration  of  the  signal  when  the  mean  time  is  zero,  and  is  a  range-doppler 
coupling  term  [6,  7]. 
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2.2  Cramer— Rao  Bound  for  a  Single  Pulse 


Similar  to  Equation  (6),  we  have  now  the  uncertainty  ellipse,  and  the  width  of  this  ellipse 
in  the  range  and  velocity  directions  is  [7] 


St 

8(f> 


1 

PVSNR 

1 

av/SNR 


(11) 

(12) 


The  aim  is  to  estimate  time-delay  and  doppler.  As  we  want  our  estimates  to  be  as  accurate 
as  possible,  we  want  the  error  measurements  to  be  as  small  as  possible  -  this  is  the  error 
variance,  and  the  minimum  error  variance  is  the  Cranrer-Rao  bound.  That  is,  our  estimate 
will  be  the  actual  value  of  the  time-delay  or  doppler,  plus/minus  an  error.  We  wish  to 
minimise  the  variance  of  this  error. 


It  is  found  that  the  minimum  error  variance 

E[(f  -  r)2] 

EM  -  a 1 02] 


1 

(13) 

f32  SNR 

1 

a2  SNR 

(14) 

with  f  and  (j)  representing  the  estimate  of  the  time-delay,  and  doppler,  respectively. 


In  the  above  the  SNR  refers  to  the  input  SNR.  Defining  SNRe  =  BT  x  (SNR)  to  be  the 
effective  output  SNR,  where  B  is  the  noise  bandwidth  and  T  the  pulse  duration,  then  we 
can  write  [8] 


E[(f  -  r)2] 

EM  -  J)?] 


1 

(32  SNRe 
1 

a2  SNRe 


(15) 

(16) 


For  a  single  pulse,  the  noise  bandwidth  is  the  inverse  of  the  sample  rate  (half  the  Nyquist 
rate) . 

In  order  to  distinguish  signals  at  different  ranges,  we  need  to  make  |x(r,  0)|  as  small  as 
possible  everywhere,  except  at  r  =  0. 

Measurement  accuracy  usually  assumes  that  one  signal  is  present,  and  depends  on  the 
signal-to-noise  ratio  and  the  behaviour  of  |y(r,  (f>)\  in  a  small  region  about  r  =  0  and 
(j)  =  (1. 

The  minimum  error  variances  given  above  require  a  knowledge  of  a  and  (3.  For  a  single 
linear  FM  pulse,  frequency  ranging  over  an  interval  A/,  with  rectangular  envelope  [7], 


a2  = 

7 T2  T 2 

3 

(17) 

/32  = 

7T2  (A/)2 

3 

(18) 

Al2  = 

tt2TA/ 

3 

(19) 
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2.3  Estimation 

We  now  consider  the  case  of  a  stationary  ESM  system.  A  waveform  has  been  emitted, 
designed  to  have  a  sharply  peaked  ambiguity  function.  As  a  result,  the  target  echo  will 
be  highly  correlated  to  the  original  signal.  The  resolution  is  determined  by  the  width  of 
the  peak. 

We  can  exploit  this  knowledge,  and  determine  angle  of  arrival,  by  using  two  passive  re¬ 
ceivers,  separated  in  space.  Knowing  the  signal  is  highly  correlated,  we  simply  correlate 
the  output  of  one  of  the  receivers  with  the  time-delayed  output  of  the  other.  This  will 
give  us  the  time-delay,  and  hence,  by  geometry,  the  angle  of  arrival. 

As  it  is  usual,  in  waveform  design,  to  work  in  the  frequency  domain,  we  do  so  here.  When 
it  is  time  to  estimate  the  time-delay,  we  will  perform  an  inverse  Fourier  transform,  and 
move  into  the  time-domain. 

The  algorithm  proceeds  as  follows. 

Step  1:  We  have  real  data  from  two  receivers,  xA{t )  from  receiver  1,  and  xB(t )  from 
receiver  2  -  two  time  series.  The  receivers  are  assumed  to  be  stationary,  and  separated  in 
space. 

Step  2:  We  obtain  the  signal  spectrum  via  a  Discrete-time  Fourier  Transform,  denoted 
XA(f)  and  XB(f). 

At  this  point  we  are  faced  with  the  problem  of  redundancy.  The  original  series  was  real, 
and  therefore  should  only  contain  positive  frequencies.  Under  a  Fourier  transform,  we 
have  generated  a  complex  function  with  positive  and  negative  frequencies.  Before  we  can 
continue,  we  must  suppress  the  negative  frequencies  [5,  6,  9]. 

It  can  be  shown  that  if  a  complex  signal  can  be  written  as  the  sum  of  a  real  part,  and  an 
imaginary  part,  and  if  these  parts  are  related  by  a  Hilbert  transform,  then  one  half  of  the 
spectrum  will  be  suppressed  [2], 

A  signal  satisfying  the  above  is  called  an  analytic  signal. 

Step  3:  Create  signals  with  no  negative  frequencies. 

Let  ip(t)  denote  the  analytic  signal. 

The  frequency  spectrum  of  ip(t)  is  [2,  5,  6] 

*(/)  =  2  X(f)  f>  0 

=  0  otherwise 

Thus,  having  obtained  the  signal  spectrum  in  Step  2,  we  simply  double  the  amplitude  for 
positive  frequencies,  and  set  all  other  amplitudes  to  zero. 

Step  4:  The  zero  doppler  Ambiguity  function  reduces  to  a  complex  cross  correlation,  in 
the  time  domain.  In  the  frequency  domain,  it  is  the  product  of  the  Fourier  transforms 
^a(/)  and 

Thus,  we  calculate  ^A(f)  \Kg(/)  -  the  Cross  Spectral  density. 
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Time  Delay 


Figure  1:  Quadratic  Interpolation 


Step  5:  We  revert  back  to  the  time  domain  via  an  inverse  Fourier  transform  of  the  power 
Spectral  density,  in  order  to  produce  the  CCF.  We  zero  pad  the  IFFT  in  order  to  inter¬ 
polate  the  CCF. 

Zero  padding  involves  interpolation  by  increasing  sample  rate.  To  increase  the  sample  rate 
by  the  factor  k,  we  need  to  calculate  k  —  1  intermediate  values  within  the  original  time 
series  [10]. 

Zero  padding  not  only  increases  the  number  of  points  in  our  new  time  series,  but  it  changes 
the  sample  rate.  If  the  original  time  series  has  N  data  points,  and  if  the  zero  padding 
factor  is  k,  then  the  new  time  series  will  have  kN  data  points.  If  the  original  time  series 
was  sampled  every  T  seconds,  then  the  new  time  series  will  be  sampled  every  T/k  seconds. 

If  we  assume  the  signal  is  a  stationary  process,  then  we  would  expect  a  single  peak  at 
the  value  of  the  time-delay.  As  a  result,  we  could  expect  the  correlation  function  to  be 
approximated  by  a  parabola  near  the  peak.  The  peak  of  this  parabola  is  found  by  parabolic 
interpolation  -  we  use  three  points,  one  on  either  side  of  the  local  maxima  (as  shown  in 
Figure  1),  to  find  the  peak,  and  estimate  the  time-delay  [11,  12,  13]. 

Using  Newton’s  ‘forward  interpolation  method’,  we  Taylor  expand  the  cross  correlation 
function  about  the  local  maxima  tm 


where 


Gw(t)  —  \Gw(tm)\  +  (t  —  tm ) — - - f  —  (t  —  tm)(t  —  tm- 1) — ^2 — 


AG0  —  \Gw(tm)\  —  \Gw(tm-i)\ 


(20) 

(21) 
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and 

A2  Go  =  |G„,(im+i)|  —  2\Gw(tm)\  +  \Gw(tm-i)\  (22) 

with  h  =  tm  —  tm- 1- 

Expanding  out  the  factors  in  the  Taylor  expansion  gives 

Gw(t)  =  c  +  t  ^ — - - — +  -^j^t2  A2Gq  (23) 

where  c  is  a  constant. 

We  see  that  if  we  write  —h  =  tm-\  —  tm,  then 

tm—  1  T  tm  =  —  (h  —  2  tm)  (24) 

and  we  can  write 

Gw(t)  =  c  +  t  ^ +  7^2^  —  2tm]A2Go^  +  2h^t2^2G°  (25) 


If  we  define 


a  = 


b  = 


1 

W2 

AG0 


/i 


A2  G0 

+  a  [h  —  2  tm] 


then 

Gw(t)  =  at2  +  bt  +  c 

with  maxima  at  the  zero  of  the  derivative  with  respect  to  the  time-delay. 
That  is,  the  estimate  of  the  time-delay,  using  quadratic  interpolation,  is 


td 


_b_ 

2  a 

i  \Gw(tm)\  —  \Gw(tm-i)\ 

\Gw(tm+i)\  —  2\Gw(tm)\  +  |GTO(im_i)| 


(26) 

(27) 

(28) 


(29) 


after  substituting  for  a  and  b. 

Substituting  this  back  into  the  above  quadratic,  gives  the  value  of  Gw  at  the  maxima, 


Go  =  Gw(td)  —  —  - — he 

4a 


(30) 


The  constant  c  can  be  written  in  terms  of  a  and  b, 

c=\Gw(tm)\ -btm-  at2m  (31) 

by  collecting  the  terms  not  shown  explicitly  in  Equation  (23)  and  using  Equation  (24). 

A  property  of  the  ambiguity  function  is  that  its  maximum  occurs  at  the  origin,  and  is 
given  by  [5] 

\x(T:4>)\  <  X(0,0)  =  2E  (32) 
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with  E  the  signals  energy.  Thus,  given  Equation  (30),  we  can  find  E. 

From  the  Taylor  expansion  of  the  ambiguity  function,  we  see  that  the  effective  bandwidth 
is 


-1  d2Gw(t) 
lx(0,0)|  dt2 


(33) 


evaluated  at  id.  That  is, 


f?  =  -|r  (34) 

and  reflects  the  curvature  of  the  ambiguity  function,  at  the  peak. 

In  terms  of  the  effective  bandwidth,  or  curvature,  we  have,  from  Equation  (29) 


id  = 


b 


(35) 


and  so,  for  fixed  Go,  the  time-delay  estimate  depends  on  the  curvature  of  the  ambiguity 
function  (at  its  peak).  The  greater  the  curvature,  the  narrower  the  ambiguity  function  at 
the  peak,  leading  to  an  improved  estimate  of  the  time-delay. 

We  see  from  the  definition  of  (3 2  given  by  Equation  (18),  that  by  increasing  the  effective 
bandwidth,  we  increase  A/,  and  as  a  result,  we  increase  chirp  rate  (if  nothing  else  changes). 

Thus  the  curvature  of  the  ambiguity  function  at  its  peak  depends  on  the  chirp  rate  -  the 
greater  the  chirp  rate,  the  narrower  the  peak. 


2.3.1  Continuous— Wave  Cramer— Rao  Bound 


The  Continuous-wave  Cramer-Rao  bound  (CW  CR-bound)  can  be  found  from  the  log 
likelihood  function  of  the  “system”  under  consideration. 

For  example,  we  consider  a  quadrature  signal  at  one  receiver,  and  its  time-delayed  form 
received  at  a  second  receiver.  Assuming  the  same  signal  amplitude  at  both  receivers,  the 
signal  functions  are 

si(in)  =  b  cos[2TTfbtn  +  7T  od2n]  +  i  b  sin[2-jrfbtn  +  vr  at2n\  (36) 

S2(tn)  =  b cos[2tt fb(tn  -  td )  +  Tra(tn  -  td )2]  +  i  b  sin[2n fb(tn  -  td)  +  ir  a(tn  -  td)2} 


where  si  refers  to  the  signal  in  receiver  1,  and  S2  refers  to  the  signal  in  receiver  2.  We  also 
define  fb  to  be  the  initial  chirp  frequency,  and  chirp 


a 


fe  "  fb 

At 


(37) 


with  fe  the  end  chirp  frequency,  and  At  the  time  it  takes  to  increase  (or  decrease)  from 
fb  to  fe. 

For  receiver  1,  define 


W(tn)  =  KM  tnj\+W(tn) 
Y\{tn)  =  tn)]  +  W(tn) 
Zi(tn)  =  Xi(tn)  +iYi(tn) 


(38) 
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with  W  ~  1V(0,(T2),  so 

Xi  (tn)  ~  N(bcos[2TTfbtn  +  natl\,a2)  =  N(/i1{n),a2) 

Yi(tn)  ~  N(b  sin[2ir  fbtn  +  Ttatn],cr2)  =  N(v\(n),a2) 

having  defined  Hi(n)  =  b  cos[2tt  fbtn  +  tt at2]  and  iq (n)  =  6sm[27r/bin  +  7raf2]. 

The  pdf  of  Z\  can  be  written 

n(r  ' |  _  r-v7T^-n  ([^i(*«)-Mi(»)]2+[yi(t«)-^i(n)]2) 

K  lj  (27ra2)N 

Similarly,  for  receiver  2,  define 

X2(fn)  =  K[S2(*n)]  +  IT(*n) 

1^2  (tn)  =  %[s2(tn)}+W(tn) 

Z2(tn)  =  X2(tn)+iY2(tn) 

so  that 

X2(tn)  ~  N(bcos[2TT fb(tn  -  td)  +  ira(tn  -  td)2],c r2)  =  N(n2(n),  a2) 
Y2{tn )  ~  fV(6.sm[27r/b(tn  -  id)  +7ra(fn  -  td)2],cr2)  =  N  (v2(n) ,  a2) 

defining  n2{n)  =  6cos[27r/fe(tn  —  td)  +  na(tn  —  td)2]  and 
u2(n)  =  b sin[2-K fb{tn  -  td)  +  n a(tn  -  td)2}. 

The  pdf  of  Z2  can  be  written 

,J7  \  ^  r-^ 7T  Y^I^([X2{tn)-^2(n)]2+\Y2{t„)-V2(n)}2  ) 

H  2j  {2tuj1)N 


(39) 


(40) 


(41) 


(42) 


(43) 


Finally,  define  the  vector 

Z(tn)  =  [Z1(tn)Z2(tn)}r  (44) 

then 

f(7' |  _  r~y^r  ^^-nL([Xi(*”)-MiW]2+[yi(*n)-^i(n)]2+[X2(tn)-M2(n)]2+[y2(tn)-^2(n)]2) 

Jl  j  (27ra2)2^ 

(45) 


The  Cramer-Rao  bound  is  obtained  from  the  Fisher  information  matrix  [14] 

-dlnf(Z)  Sin f(Z)- 


Jij  —  E 


da ; 


da^ 


(46) 


with  1  <  i,j  <  k,  for  k  parameters.  In  the  CW  case,  the  cq  are  elements  of  the  set 
(./),,  b,  t,  a),  so  k  =  4.  The  symbol  E  in  the  above  Equation  represents  the  expectation. 


Performing  the  differentiation,  and  expanding,  we  find  that 


J, 


ij 


1 


o“ 


N- 1 


E 

n= 0 


+rx  ( 


+rY  ( 


dm(n)  dni(n) 


l  dXi  dxj 
d/jti(n)  d/j2(n) 
dXi  dxj 
dv\(n)  dv2{n) 
dXi  dxj 


+ 

+ 


d/J,2(n)  dn2(n) 

*  r\ 


+  dxi 

dXj 

dn2(n) 

dm  (n) 

dxi 

dXj  J 

dis2(n)  dv\ (n)J 

dXi 

dXj  J\ 

dv\(n)  dv\(n)  du2(n)  dv2{ri) 
+  dXi  dxj  +  dxi  dxj 


l<i,j<4  (47) 
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for  the  Fisher  information  matrix,  with  xi  =  fb ,  X2  =  b,  X3  =  td  and  X4  =  a- 

A  consequence  of  the  signal  model  we  have  chosen  is  that  there  will  be  cross  terms  in  the 
Fisher  information  matrix.  In  the  above,  rx  and  ry  are  the  correlation  coefficients  for  X 
and  Y .  That  is, 

coviX  i,X2) 

T  Y  —  —  — 

\/var(Xi)var{X2) 

with  —  1  <  rx  <  1-  Similarly  for  ry. 

In  what  is  to  follow,  rather  than  calculate  the  correlation  coefficients,  we  consider  two 
special  cases:  rx  =  ry  =  1  for  high  SNR,  when  we  can  expect  some  correlation,  and 
rx  =  ry  =  0  for  low  SNR,  and  little  correlation.  We  justify  this  by  considering  the 
number  of  data  samples  receiver  2  lags  behind  receiver  1.  For  small  time-delays,  up  to 
20  ns,  say,  receiver  2  will  only  be  a  few  samples  behind  receiver  1  (eg,  sampling  at  5  ns). 
For  high  SNR,  we  can  expect  the  waveform  in  receiver  2  to  follow  that  in  receiver  1,  except 
near  the  peaks  and  troughs  for  sinusoids,  and  the  beginning  and  end  of  the  frequency  ramp 
for  chirped  signals.  For  low  SNR,  we  can  not  expect  this. 

A  consequence  of  the  non-linear  time  terms  in  the  signal  model  is  that  the  calculation  of 
this  matrix,  and  the  inverse,  becomes  quite  complicated.  For  the  case  where  there  is  no 
chirp  (the  three  dimensional  case),  we  obtain  expressions  for  the  variance  of  /  and  b  which 
are  similar  (up  to  a  factor)  to  those  found  by  Rife  and  Boornstyn  [14],  The  variances  are 
smaller  by  a  factor  of  4,  due  to  the  fact  we  are  considering  the  two  receiver  case.  The 
variance  of  td  shows  it  depends  explicitly  on  /  (again  we  are  only  considering  the  two  cases 
of  rx  =  ry  =  1  and  rx  =  ry  =  0). 

The  four  dimensional  case  shows  that  the  variances  of  the  parameters  under  consideration 
are  dependent  on  the  estimates  of  those  parameters. 

For  example,  at  high  SNR  (rx  =  ry  =  1), 

/o  l\  2  N—l 

J33  =  (  - —  )  J2ifb  +  a(tn  -  id)2}  (49) 

which  we  see  is  a  function  of  td,  the  estimate  of  the  time-delay. 

For  the  case  where  fb,  b  and  a  are  known,  the  Cramer-Rao  bound  for  the  variance  of  the 
time-delay  would  simply  be  the  inverse  of  J33.  For  the  case  where  they  are  all  unknown, 
the  variance  will  be  a  function  of  J33. 

We  should  point  out  that  this  dependence  is  not  entirely  unexpected.  Rife  and  Boornstyn 
derive  expressions  for  the  variance  of  frequency  and  phase,  both  of  which  show  amplitude 
dependence  [14]. 
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Source  i 


Antenna  1 


Antenna  2 


Cross  correlator 


Figure  2:  Geometry  of  receiver  system 


2.3.2  Simulations 


2.3.2. 1  Continuous  Wave 

The  algorithm  described  above  has  been  coded  in  Matlab.  A  continuous  wave  (CW)  signal 
is  linearly  swept  from  fb  =  5  MHz  to  fe  =  35  MHz,  over  At  =  1  milli  second,  to  produce  a 
sawtooth  waveform  over  the  sample  time.  This  sawtooth  waveform  will  allow  us  to  observe 
the  effects  of  frequency  wrap  around  on  our  estimates. 

Figure  2  shows  the  geometry  of  the  receiver.  From  this  we  have  the  time-delay 


td  = 


L 

c 


sin  9 


(50) 


where  c  is  the  speed  of  light,  and  9  is  the  angle  of  arrival. 


We  show  the  instantaneous  frequency  in  Figure  3,  plotted  using  the  Matlab  specgram 
function  [15]. 


The  instantaneous  frequency  is  proportional  to  the  time  derivative  of  the  phase  of  the 
signal.  Namely, 


1  d  <j) 
2n  dt 


(51) 


with 

cj>  =  2tt  (ft  +  at2 /2) 


(52) 
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Time  micro  seconds 

Figure  3:  Instantaneous  frequency  of  chirp  waveform  over  sample  time 

An  example  of  the  cross  spectral  density  is  shown  in  Figure  4.  It  shows  clearly  the  effect 
of  noise  and  the  sawtooth  shape  of  the  instantaneous  frequency.  In  this  Figure,  SNR  is 


10  dB. 
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Figure  4:  Cross  Spectral  Density  -  CW 
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In  Figures  5  to  16  we  plot  histograms  of  the  estimate  of  the  time-delay.  The  simulation 
used  to  estimate  the  time-delay  ran  over  1024  iterations,  with  218  data  points,  and  a 
sample  rate  of  5  ns.  This  number  of  points  is  used  for  simulation  purposes  only. 

Although  this  report  is  mainly  concerned  with  time-delay  estimation,  angle  estimation 
will  be  required  when  we  come  to  the  EWRD  Durex  trial  data.  As  such,  we  also  give 
angle  estimates  in  the  following  figures.  Angle  estimates  are  obtained  from  the  time-delay 
estimate  after  each  iteration,  by  inverting  Equation  (50),  and  replacing  tfj  with  ta¬ 
in  Figures  5  to  10,  we  fix  the  antenna  separation  at  L  =  10  m,  and  vary  0  and  SNR.  In 
Figures  11  to  16,  we  fix  the  antenna  separation  at  L  =  30  m,  and  vary  9  and  SNR.  We  use 
different  true  time-delays  for  both  separations  in  order  to  consider  a  wide  range  of  values. 

The  estimate  of  the  time-delay  was  taken  to  be  the  value  with  the  maximum  count  number 
in  the  histogram.  This  value  was  then  used  in  the  estimate  of  the  angle  of  arrival. 
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Figure  5:  Time-delay  estimate.  L  =  10  m,  td  =  5.8  n  sec,  CW 
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Figure  6:  Angle  estimate.  L  =  10  m,  9  =  10°,  CW 
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We  see  bias  in  Figures  5  to  10  for  antenna  separation  of  10  nr.  We  also  note  the  large 
peaks  at  SNR  below  -10  dB.  Although  we  may  be  tempted  to  use  these  values,  we  must 
take  into  account  the  large  variance. 

Again  we  note  that  these  figures  were  generated  from  a  simulation  of  1024  iterations  with 
218  data  points.  As  such  these  estimates  will  be  more  accurate,  and  the  variances  smaller 
than  we  would  expect  from  real  data. 
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Figure  1 7;  Variance  versus  SNR.  L  =  10  m,  0  =  30° 
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*  Variance  of  time-delay  estimate 
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Figure  18:  Variance  versus  SNR.  L  =  10  m,  6  =  10° ,  CW 


The  above  simulations  suggest  a  bias  in  the  estimate  of  time-delay,  perhaps  as  a  result 
of  the  frequency  wrap  around.  For  a  receiver  separation  of  10  nr,  we  see  the  difference 
between  the  actual  time-delay  and  the  estimated  time-delay  increasing  as  the  angle  of 
arrival  increases.  For  an  angle  of  arrival  of  80°,  and  an  SNR  of  10  dB,  we  see  the  estimated 
time-delay  differs  from  the  actual  time-delay  by  approximately  5  ns,  or  25°.  The  actual 
time-delay  in  this  case  being  32.8  ns. 

For  a  receiver  separation  of  30  nr,  the  error  is  small  and  appears  constant.  This  suggests 
small  antenna  separations  may  not  be  appropriate  for  CW  signals. 
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2. 3. 2. 2  Single  Pulse 

In  Figures  21  to  32  we  plot  histograms  of  the  estimate  of  the  time-delay.  The  simulation 
used  to  estimate  the  time-delay  ran  over  1024  iterations,  with  218  data  points.  Sample 
rate  is  5  ns.  As  with  the  CW  case,  we  also  give  angle  estimates.  Angle  estimates  are 
obtained  from  the  time-delay  estimate  after  each  iteration,  by  inverting  Equation  (50), 
and  replacing  td  with  td- 

In  Figures  21  to  26,  we  fix  the  antenna  separation  at  L  =  10  nr,  and  vary  9  and  SNR.  In 
Figures  27  to  32,  we  fix  the  antenna  separation  at  L  =  30  m,  and  vary  9  and  SNR. 

The  estimate  of  the  time-delay  was  taken  to  be  the  value  with  the  maximum  count  number 
in  the  histogram.  This  value  was  then  used  in  the  estimate  of  angle. 

Figure  19  was  obtained  using  the  Matlab  specgram  function  [15] .  There  is  no  wrap-around. 
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Figure  1 9:  Instantaneous  frequency  of  chirp  waveform  over  sample  time  -  single  pulse 
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Figure  20:  Cross  Spectral  Density  -  single  pulse 
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In  contrast  to  the  CW  case,  the  single  pulse  simulations  give  estimates  which  are  very  close 
to  the  actual  time-delay  and  angle  of  arrival  values  down  to  0  dB.  However  the  spread 
of  values  increases  significantly  as  SNR  decreases.  We  also  assume  that  the  correct  signal 
has  been  detected  at  such  a  low  SNR. 

Again  we  note  these  estimates  are  expected  to  be  more  accurate,  and  the  variances  smaller 
than  we  would  expect  from  real  data. 


32 


DSTO-TR  1705 


3  EWRD  Durex  Trial  Data 

We  now  apply  the  techniques  described  earlier  to  data  obtained  from  the  EWRD  Durex 
trial  of  17  July,  2002  [1], 

In  Figure  33,  we  see  the  true  time-delay  as  a  function  of  angle.  This  figure  shows  that 
for  large  angle,  small  changes  in  time-delay  represent  large  changes  in  angle.  For  this 
reason,  good  estimates  will  be  found  where  the  curve  is  not  near  the  regions  of  maximum 
curvature  -  that  is,  for  —60°  <  9  <  60°. 

The  samples  were  of  pseudo-random  noise  with  20  MHz  bandwidth. 

The  receiver  was  dual  channel,  with  8  bit  ADC,  sampling  at  200  MSps  at  50  MHz  IF. 
The  two  receivers  were  29.5  m  apart,  with  the  straight  line  between  antennas  being  0° 
magnetic  North. 

The  two  receivers  were  not  calibrated. 

The  transmitter  was  moved  parallel  to  the  receivers,  with  angle  being  measured  using  a 
compass.  The  last  series  of  data  were  taken  due  south  of  the  receivers. 

The  geometry  of  this  setup  is  shown  in  Figure  34. 

We  compare  the  geometry  of  the  Durex  trial  with  that  given  by  Figure  2,  and  observe  that 
a  mapping  between  the  two  conventions  is  needed.  We  see  that  the  Durex  trial  definition 
of  90°,  is  the  theoretical  0°.  Similarly,  135°  becomes  —45°,  150°  becomes  —60°,  and  22° 
becomes  68°. 

We  note  that  negative  angles  refer  to  receiver  1  lagging  receiver  2. 


Figure  33:  Time-delay  as  a  function  of  angle 
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Figure  34:  Receiver  -  Transmitter  geometry,  Durex  trial 
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In  Table  3.1,  we  present  trial  estimates  of  the  angle.  The  angle  given  in  column  4  is  that 
found  by  a  compass  reading.  Table  3.2,  gives  estimates  of  the  time-delay.  The  time-delay 
given  is  the  expected  delay  given  the  angle  and  receiver  separation  of  29.5  m. 

Figure  35  shows  time-delay  as  a  function  of  angle.  We  have  also  plotted  estimates  of  the 
time-delay,  given  in  Tables  3.1  and  3.2. 

Unfortunately,  due  to  the  nature  of  the  data  collection,  we  are  unable  to  comment  on  the 
accuracy  of  the  time-delay  estimation  method  using  the  Durex  trial  data. 

Table  3.1:  Angle  estimates  from  Durex  trial  data 


File  name 

Tx  Power 
(dBm) 

Rx  attenuation 
(dB) 

Angle  (°) 

Estimated  Angle  (°) 
(averaged) 

ADC001  -  ADC009 

39 

25 

0 

-4.7 

ADC010  -  ADC019 

39 

25 

-45 

-35.2 

ADC020  -  ADC029 

0 

0 

-45 

-31.2 

ADC030  -  ADC039 

0 

0 

-60 

-41.8 

ADC040  -  ADC049 

39 

25 

-60 

-44.8 

ADC080  -  ADC089 

39 

25 

0 

2.2 

ADC090  -  ADC099 

39 

25 

45 

48.8 

ADC100  -  ADC109 

39 

25 

68 

77.3 

ADC110  -  ADC119 

39 

25 

-90 

-62.4 

Table  3.2:  Time-delay  estimates  from  Durex  trial  data 


File  name 

Tx  Power 
(dBm) 

Rx  attenuation 
(dB) 

Delay  (ns) 

Estimated  Delay  (ns) 
(averaged) 

ADC001  -  ADC009 

39 

25 

0 

-8.1 

ADC010  -  ADC019 

39 

25 

-69.5 

-56.7 

ADC020  -  ADC029 

0 

0 

-69.5 

-50.9 

ADC030  -  ADC039 

0 

0 

-85.2 

-65.5 

ADC040  -  ADC049 

39 

25 

-85.2 

-69.3 

ADC080  -  ADC089 

39 

25 

0 

3.8 

ADC090  -  ADC099 

39 

25 

69.5 

74 

ADC100  -  ADC109 

39 

25 

91.2 

95.9 

ADC110  -  ADC119 

39 

25 

-98.3 

-87.1 
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Figure  35:  Durex  trial,  time-delay  versus  angle. 
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4  Time— Domain  Filtered  Cross  Spectral  Density 

In  the  previous  chapter,  the  algorithm  for  estimating  time-delay  was  :  start  in  the  time 
domain  with  data  from  two  spatially  separated  receivers.  We  move  into  the  frequency 
domain  by  performing  a  Fourier  transform  on  each  data  set.  We  suppress  negative  fre¬ 
quencies  by  generating  an  analytic  signal  via  a  Hilbert  transform,  and  we  multiply  their 
Fourier  transforms  to  produce  the  cross  spectral  density. 

We  then  apply  an  inverse  Fourier  transform  to  this  density  to  move  back  into  the  time 
domain,  to  give  the  cross  correlation  function.  We  time  filter  this  function,  to  remove 
unnecessary  data  -  away  from  the  peak  -  and  perform  a  quadratic  interpolation  on  the 
filtered  function  to  give  us  the  time-delay  estimate. 

In  [16],  the  cross  correlation  is  performed  in  the  time  domain,  and  the  cross  spectral 
density  is  obtained  via  Fourier  transform.  The  time-delay  is  estimated  in  the  frequency 
domain  from  the  phase  slope  of  the  cross  spectral  density. 

Then  a  spread-spectrum  signal  arriving  at  two  spatially  separated  receivers  is  cross  corre¬ 
lated.  The  cross-correlation  output  will  be  the  sum  of  1)  a  time-delayed  auto-correlation 
of  the  signal,  2)  the  cross  correlation  of  the  signal  with  noise,  and  3)  the  cross-correlation 
of  the  noise. 

The  noise  cross-correlation  component  will  be  spread  uniformly  over  the  entire  cross¬ 
correlation  function.  The  auto-correlation  component  will  peak  near  the  center  of  this 
function.  The  shift  from  the  center  represents  the  time-delay.  However,  noise  will  ensure 
the  peak  does  not  occur  at  the  actual  time-delay,  so  it  is  not  possible  to  just  “read  off” 
the  value  from  the  plot,  with  reasonable  accuracy. 

Knowing  that  the  auto-correlation  component  will  be  concentrated  over  a  region  tens  of 
nanoseconds  wide,  we  can  apply  a  low  pass  time-domain  filter.  An  FFT  applied  to  this 
filtered  function  leads  to  the  time-domain  filtered  cross  spectral  density.  It  also  leads 
to  phase  wrapping.  The  phase  slope  of  this  density,  across  the  signal  bandwidth,  is  the 
estimate  of  the  angle  of  arrival  [16]. 

If  the  signal  is  shifted  by  an  amount  td,  then  the  rate  of  change  of  phase  with  frequency  is 
proportional  to  td-  Once  the  time-domain  filtered  cross  spectral  density  has  been  obtained, 
the  phase  slope  can  be  estimated  using  linear  regression.  From  this,  we  estimate  time- 
delay. 

We  mentioned  that  the  application  of  the  FFT  leads  to  phase  wrapping.  Therefore,  we 
need  to  unwrap  the  phase  before  we  make  a  straight-line  fit.  This  requires  any  phase 
errors  at  each  sample  to  be  small.  It  has  been  shown  that  the  phase  errors  due  to  noise 
are  small  if  SNR  >  13  dB  [16,  17]. 

The  angle  of  arrival  estimation  error  is  found  to  be  similar  in  form  to  Equation  (A14).  It 
was  found  that  high  accuracy  can  be  obtained  (fractions  of  a  degree),  with  large  errors 
found  for  \0\  >  60°.  The  estimate  is  also  shown  to  be  independent  of  the  receiver  noise 
bandwidth  and  the  “size”  of  the  time-domain  filtering,  providing  it  is  large  enough  to 
enable  detection. 
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5  Discussion 

The  preceding  was  a  technical  report  on  the  usefulness  of  exploiting  the  waveform  design 
of  highly  correlated  LPI  signals,  in  order  to  estimate  time-delay,  and  as  a  result,  angle  of 
arrival.  The  model  consisted  of  two  stationary  receivers  spaced  several  metres  apart.  The 
signal  source  is  also  assumed  to  be  stationary. 

For  the  simulations,  we  considered  a  chirped  signal.  The  time  taken  from  the  lowest  to 
highest  frequency  value  was  1  ms.  The  signal  was  sampled  evenly,  the  sample  rate  being 
5  ns,  at  various  SNRs,  and  the  time-delay  estimated,  as  described  in  the  report.  Quadratic 
interpolation  was  used  to  obtain  sub  sample  rate  estimates.  For  the  purpose  of  simulation, 
we  used  218  data  points,  and  1024  iterations.  As  a  result,  we  can  expect  our  estimates  to 
be  more  accurate,  and  our  variances  smaller  than  would  be  expected  from  real  data. 

We  considered  the  case  of  a  continuous  wave  (CW)  signal,  and  a  single  pulse.  The  CW 
case  allowed  us  to  observe  the  effects  of  frequency  wrap  around  on  the  estimates.  For  the 
CW  signal,  the  sample  time  is  just  over  1  ms,  allowing  for  frequency  wrap  around.  For 
the  pulsed  signal,  sample  time  is  the  pulse  duration,  which  is  1  ms. 

We  started  the  simulation  with  a  CW  signal,  a  receiver  separation  of  10  m,  and  a  true 
time-delay  of  5.8  ns,  representing  a  true  angle  of  arrival  of  10°.  As  we  found  in  Section 
2.3.2,  our  time-delay  estimates  were  off  by  approximately  1  ns,  or  approximately  2°.  That 
is,  at  an  SNR  of  10  dB,  and  a  true  time-delay  of  5.8  ns  (10°),  we  estimated  the  time-delay 
to  be  4.8  ns,  with  a  variance  of  0.02  ns.  This  time-delay  estimate  represents  an  angle  of 
arrival  estimate  of  8.2°. 

Even  at  -10  dB  SNR,  we  estimated  the  time-delay  to  be  4.9  ns.  However,  the  variance 
was  1.5  ns,  with  the  spread  of  estimates  far  greater  than  those  found  at  10  dB. 

We  then  increased  the  true  time-delay  (or  angle  of  arrival)  for  the  same  antenna  separation 
of  10  nr.  We  found  the  estimates  differed  from  the  true  values  by  an  increasing  amount. 
That  is,  at  10  dB,  for  a  true  time-delay  of  16.7  ns  (30°),  we  estimated  13.7  ns  (24.1°).  For 
a  true  time-delay  of  32.8  ns  (80°),  we  estimated  27.5  ns  (55.7°). 

We  then  increased  antenna  separation  to  30  m,  and  repeated  the  simulations  but  with 
different  time-delays.  In  this  case,  at  10  dB,  we  found  the  difference  between  the  true 
time-delay  and  the  estimated  time-delay  to  be  approximately  2  to  3  ns. 

At  10  dB,  for  a  true  time-delay  of  17.4  ns  (10°),  we  estimated  14.3  ns  (8.2°).  For  a  true 
time-delay  of  70.7  ns  (45°),  we  estimated  73.8  ns  (47.6°).  For  a  true  time-delay  of  94  ns 
(70°),  we  estimated  95.8  ns  (73.4°). 

As  SNR  decreased  we  found  the  estimates  to  be  constant,  although  the  spread  of  estimates 
increased. 

The  above  suggests,  for  CW  signals,  there  may  be  an  optimum  antenna  separation  for 
which  estimates  are  reasonably  accurate,  and  errors  are  small  for  all  angles  of  arrival. 

Repeating  the  above  for  a  pulsed  signal  produced  results  more  accurate  than  the  above. 

For  a  receiver  separation  of  10  m,  and  a  true  angle  of  arrival  of  30°,  representing  a  true 
time-delay  of  16.7  ns,  we  estimated  the  time-delay  to  be  16.7  ns,  with  a  variance  of  0.03 
ns.  This  time-delay  estimate  represents  an  angle  of  arrival  estimate  of  30°. 
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For  a  receiver  separation  of  10  m,  and  a  true  angle  of  arrival  of  10°,  representing  a  true 
time-delay  of  5.8  ns,  we  estimated  the  time-delay  to  be  5.8  ns,  with  a  variance  of  0.03  ns. 
This  time-delay  estimate  represents  an  angle  of  arrival  estimate  of  10°. 

This  analysis  suggests  that,  even  with  SNR  at  -10  dB,  this  method  gives  a  reasonable 
measure  of  the  time-delay,  or  angle  of  arrival.  However,  the  spread  of  estimates  may 
negate  the  usefulness  of  such  estimates,  as  well  as  the  need  to  actually  observe  a  signal  at 
such  a  low  SNR. 

We  have  also  tested  this  approach  against  pseudo-random  noise  generated  during  the 
EWRD  Durex  trial.  Due  to  the  data  obtained,  we  were  unable  to  comment  on  the  accuracy 
of  our  estimation  method. 

Possible  future  research  : 

When  estimating  the  time-delay,  we  made  use  of  a  quadratic  interpolation,  in  order  to 
achieve  sub  sample  rate  accuracy.  This,  however,  leads  to  a  biased  estimate  [11]. 

The  pulsed  case  seems  to  suggest  this  bias  will  be  small.  This  may  not  be  so  for  the  CW 
case,  although  the  frequency  wrap  around  effect  may  well  dominated  any  bias.  This  could 
be  investigated  further. 

It  is  also  possible  that  this  bias  can  be  reduced  by  ‘windowing’. 

The  approach  described  in  this  report  should  be  tested  against  Pilot-like  signals  [22],  with 
accurate  ground  truth.  This  will  also  suggest  the  amount  of  data  needed  for  accurate 
estimation  in  “real  world”  situations. 

Our  approach  to  calculating  the  Cramer-Rao  bounds  was  a  “brute  force”  approach,  with 
a  number  of  assumptions. 

Is  there  a  model  in  which  the  signal  functions  in  the  two  receivers  are  orthogonal?  This 
would  remove  the  cross  terms  from  the  log  likelihood  function. 

The  obvious  case  would  be  to  have  the  signal  in  the  second  receiver  90°  out  of  phase  with 
the  signal  in  the  first  receiver.  This,  however,  is  highly  unlikely  to  occur. 

Re-writing  the  problem  in  terms  of  analytic  functions  may  be  useful  [2]. 

One  obvious  generalisation  to  the  theory  given  in  this  report  would  be  to  include  a  second 
signal  source,  and  consider  the  detection  and  estimation  problem  [23,  24]. 

Another  would  be  to  consider  moving  sources  and/or  receivers. 
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Appendix  A  Optimum  Passive  Bearing 

Estimation 

The  paper  of  MacDonald  and  Schultheiss  [18]  considers  the  problem  of  ‘optimum  passive 
bearing  estimation  by  means  of  a  passive  sonar  array’.  They  consider  M  hydrophones 
arranged  in  a  linear  array.  The  following  uses  M  =  2. 

Four  assumptions  are  made 

1.  The  signal  and  noise  are  stationary  Gaussian  processes. 

2.  The  signal  source  is  at  such  a  distance  from  the  receivers,  that  the  received  wavefront 
is  planar  over  the  dimensions  of  the  receiving  system. 

3.  Noise  is  statistically  independent  from  each  other  and  the  signal. 

4.  The  observation  time  T  is  much  larger  than  the  signal  and  noise  correlation  time 
and  the  travel  time  of  the  wavefront  across  the  receiving  system. 

They  start  by  finding  the  Cramer-Rao  lower  bound  of  the  unbiased  estimators  mis  error. 
They  them  compare  this  to  the  rms  error  from  a  split -beam  tracker.  The  comparison 
suggests  that  the  split-beam  tracker  is  close  to  optimal. 


A.l  Cramer-Rao  Bound 


Over  the  observation  time  T,  we  have  (continuous  time)  data  received  at  antenna  A 


XA(t)  =  SA(t)  +  nA(t) 


(Al) 


and  at  antenna  B 

xB(t)  =  sA(t  -  td)  +  nB{t)  (A2) 

We  see  explicitly  that  the  data  received  at  antenna  B  is  the  time-delayed  data  received  at 
antenna  A.  The  aim  is  to  estimate  the  time-delay  td,  denoted  td,  and  the  variance  of  this 
estimate,  denoted  aj  ,  given  the  data. 

td 

The  minimum  variance  of  an  unbiased  estimate  is  the  Cramer-Rao  bound,  and  is 


atd> 


-1 


E 


d2L(x) 
dt;2d  J 


(A3) 


where  E  denotes  the  expectation,  and  L(x )  is  the  likelihood  function,  or  more  commonly, 
the  log  likelihood  function.  The  likelihood  function  may  be  found  via  the  Fourier  coeffi¬ 
cients  of  the  received  data. 


We  are  able  to  perform  Fourier  transforms  on  this  data  in  order  to  create  a  new  vector. 
Let 


XA(k)  = 

J  xA{t)  e~i2*kT  dt 

XB(k)  = 

J  xB(t)e~i2nk^  dt 

k  =  l...N 

(A4) 
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Then  a  new  2 N  x  1  vector  can  be  written 


X  =  [XA(1)XA(2) . . .  XA(N)XB(  1)Xb(2)  . . .  Xb(N)}'  (A5) 


where  the  ’  superscript  refers  to  the  transpose  operation. 

The  complex  probability  density  function  (pdf)  for  X  can  be  written 


p(X)  = 


-Eli  U)*aU)+x *b  (j)xB(j)) 


ttn  7TjJ=iaj 


-X*  T.~lX 


7T 


2N  lS| 


(A6) 


where  the  *  refers  to  complex  conjugation.  Here  |X|  is  the  determinant  of  the  2 N  x  2N 
covariance  matrix  X  (not  4iV  x  4iV  as  X(k)  =  X*(T  —  k),  so  not  all  the  vectors  are 
independent). 

As  [18]  dealt  with  bearing  estimation,  we  continue  the  analysis  in  terms  of  the  angle  of 
arrival  8. 

The  covariance  matrix  X  =  E[X*X'],  and  the  matrix  elements  are 

rT/ 2  rT/2 


E[X;(wi)XJun)]  = 


l-T/2  J—T/2 


el(u)tt  dtdU 


(A7) 


where  x(t)  =  x(t,8)  is  a  function  of  the  angle  to  be  estimated.  That  is,  xp(t,6 )  = 
sA{t  —  tdp)  +  np(t),p  =  A,B,  with  tdA  =  0  and  tdB  =  (d/c)smO.  Here  d  is  the  antenna 
separation  and  c  the  speed  of  light. 

At  this  point  we  substitute  for  x(t)  in  the  expectation,  and  write  the  derived  auto¬ 
correlation  functions  for  the  signal  and  noise  in  terms  of  the  corresponding  spectral  func¬ 
tions,  remembering  that  the  signal  and  noise  are  assumed  to  be  stationary  processes. 

Working  through  the  analysis,  it  is  found  that 

E[x;(ut)x'q(un)]  =  T[S{un)+N{un)5pq]e^dq-dp^  l  =  n 

=  0  l^n  (A8) 


where  S(lo)  and  N(lo)  are  the  signal  and  noise  spectral  functions,  respectively,  and  5pq  is 
the  Kroneker  delta  function. 

/ 

The  result  of  this  is  to  separate  out  the  frequency  components  in  the  quadratic  X*  X-1A, 
and  reduce  the  original  matrix  with  2 N  x  2N  elements  to  a  matrix  with  4A"  non-zero 
elements. 

As  a  result  of  this  reduction,  we  can  reduce  the  complexity  of  the  matrix  multiplications 

/  / 

by  writing  the  quadratic  X*  X_1A,  as  the  sum  of  N  matrix  products  of  the  form  A*  BA, 
where  A  is  1  x  N,  and  B  is  N  x  N.  Each  matrix  product  involves  a  particular  frequency. 

As  a  simple  example,  suppose  we  have  M  =  2  receivers  and  N  =  2  frequencies.  Then 
our  data  vector  would  be  X'  =  [Xi(u>\)Xi(uj2)X2(uji)X2(uj2)\-  The  general  NM  x  NM 
covariance  matrix  would  look  like 
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£n(H)  £n(12)  £12(11)  £12(12) 

£n(21)  £n(22)  £12(21)  £12(22) 

£21(11)  £21(12)  £22(11)  £22(12) 

£21(21)  £21(22)  £22(21)  £22(22) 

where  the  subscripts  refer  to  vectors  X\  and  X2,  and  the  values  in  the  parenthesis  refer 

to  the  frequencies  u\  and  L02 

From  the  above  result,  only  those  frequencies  which  are  equal  have  non-zero  values,  so 

this  matrix  reduces  to  the  form 

'  £ii(ll)  0  £12(11)  0 

0  £n(22)  0  £12(22) 

£21(11)  0  £22(11)  0 

0  £21(22)  0  £22(22)  _ 

The  inverse  will  have  a  similar  form,  namely 
<722(11)  0  -<7i2(ll) 

^-1  _  1  0  022(22)  0 

*  £721  (II)  0  £7ll(ll) 

0  -£721  (22)  0 


It  is  now  possible  to  find  the  Cramer-Rao  bound  on  the  variance  of  9.  In  expanding  out 
the  matrix  products  in  the  above  quadratic,  we  find  that  there  are  terms  that  do  not 
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depend  on  the  time-delay,  or  angle  of  arrival.  We  find  that  we  can  write  the  log  likelihood 
function  as 

rp  N 

l (9)  =  —  ]T  S(k)  [X*(k)X2{k)e-ik(dl-d^  +  X*(k)X1(k)ei^dl~d2'>}  (A9) 

k= 1 

and  therefore 

-  d2)  cose  y  k  S(k)  [-X*(k)X2(k)e-ik^dl-d2'>  +  X^^X^e^1-^} 

09  CISI  /,■  i 

(A10) 

It  is  straightforward  to  show  that  the  expectation  of  this  derivative  is  zero. 

Noticing  that  terms  in  the  second  derivative  vanish  when  we  take  the  expectation,  we  find 
that  the  only  term  necessary  from  the  second  derivative  is 


82l{6) 

86 2 


T 

^2|— (di-d2)2  COS2  0 


and  the  expectation  of  this  is 


N 


E 


k2  s(k)  [xi^x^ky-^-^+x^x^ky^-y 

(All) 


E 


82i(oy 

86'2 


2  T2 


(di 


N 

d2)2  cos2  6  y  k2  S2(k) 

k=l 


(A12) 


We  also  find  that 

T2  S2(k)  _  S2(k)/N2(k) 
|Sj  “  1  +  25(jfe)/JV(fc) 


(A13) 


From  all  of  this,  we  find  the  Cramer-Rao  bound  on  the  variance  of  the  estimate  of  the 
angle  of  arrival 

N 

—  d2)2  cos2  9  'y  k2 
k= 1 


;(di 


S2(k)/N2(k) 


n-i 


l  +  2S{k)/N(k)\ 


(A14) 


The  general  form  is  given  by  Equation  (16)  in  [18]. 

This  is  the  lower  bound  of  the  mis  error  of  an  unbiased  estimator.  If  the  actual  mis  error 
achieves  this  bound,  we  have  an  optimal  estimator.  This  will  not  always  be  the  case  in  the 
real  world,  in  which  case  we  must  deal  with  ‘close  to  optimal’  estimators,  or  sub-optimal 
estimators. 
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Appendix  B  Generalised  Cross  Correlation 

Method 

This  is  a  Maximum  Likelihood  estimator  of  the  time-delay  between  signals  received  at 
two  spatially  separated  receivers  [19,  20,  21]. 

The  estimator  is  a  pair  of  pre-filters  followed  by  a  cross  correlator.  The  maximum  of  this 
correlator  in  the  time  domain,  is  the  estimate  of  the  time-delay. 

If  we  have  a  signal  from  an  unknown  source,  arriving  at  two  receivers  separated  in  space, 
then  the  received  data  can  be  written 

xi  (t)  =  si(t)+ni(t) 

x2(t)  =  a  si(t  -  td)  +  n2(t)  (Bl) 

where  tj  is  the  time-delay,  signal  s\(t)  is  uncorrelated  with  the  noise  roi(t)  and  n2(t),  and 
a  is  attenuation. 

A  common  method  for  estimating  time-delay  is  to  compute  the  cross  correlation  function 
(CCF), 

RXiX2(t)  =  E[xi(t)  x2(t  -  t)]  (B2) 

where  E  represents  the  expectation.  The  r  which  maximises  this  function  is  the  estimate 
of  the  time-delay. 

As  we  will  only  have  access  to  a  finite  duration  signal,  the  LHS  of  Equation  (B2)  is  also 
an  estimate. 

Via  a  Fourier  Transform,  we  are  able  to  write  the  CCF,  Rx  ia,2(r),  iR  terms  of  the  cross 
power  spectral  density  function  (PSD)  GXlX2(f), 

/OO 

GXlX2(f)ei2^T  df  (B3) 

-OO 

At  this  point  we  introduce  the  filters,  and  remember  that  we  only  have  an  estimate  of  the 
CCF  and  PSD.  In  terms  of  the  filtered  data  y±  and  y2,  we  have 

/OO 

H\{f)  H2(f)  GXlX2(f)ei2nfT  df  (B4) 

-OO 

with  ip(f)  =  H\  (/)  iL|(/)  being  defined  as  a  generalised  frequency  weighting.  The  weight¬ 
ing  may  be  selected  to  optimise  certain  criteria  -  for  example,  to  pass  signals  for 
which  the  signal-to-noise  ratio  (SNR)  is  high 

We  use  Equation  (B4)  to  evaluate  time-delay  estimates. 

Given  our  signal  model  from  Equation  (Bl),  we  can  write  down  the  CCF  (in  the  time 
domain).  A  Fourier  Transform  gives  us  the  PSD  (in  the  frequency  domain).  This  is  found 
to  be  the  sum  of  the  signal  PSD,  GSlSl(f),  multiplied  by  a  complex  exponential,  and  the 
noise  PSD,  Gnin2(f).  For  uncorrelated  noise,  Gnin2(f)  =  0. 

A  multiplication  in  the  frequency  domain  may  be  written  as  a  convolution  in  the  time 
domain,  and  assuming  uncorrelated  noise,  we  have 

H1112W  —  (%RXix2  (t")  *  d(t  td)  (B5) 
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The  effect  of  this  convolution  is  to  spread  out  the  peak  of  the  CCF.  A  problem  with  this 
is  when  we  have  multiple  time-delays.  This  spreading  effect  may  cause  peaks  to  overlap. 

Therefore,  the  frequency  weighting  ip(f)  is  chosen  to  ensure  large  sharp  peaks  in  the  CCF. 
Unfortunately  this  leads  to  estimates  which  are  sensitive  to  (finite  time)  errors,  particularly 
at  low  SNR. 

Let  Xi(u)),i  =  1,2  denote  the  Fourier  Transform  of  the  received  data  Xi(t). 

If  the  duration  of  the  signal,  T,  is  large  compared  to  the  absolute  value  of  the  time-delay, 
\td\,  and  the  correlation  time  of  the  signal  correlation  function,  then  we  are  able  to  make 
use  of  Equation  (A8),  and  uncouple  the  frequency  components. 

Following  through,  defining  the  vector 

X(k)  =  [A1(fe)X2(fe)]'  (B6) 

we  have  the  conditional  pdf 

p(X\Q)  =  (B7) 


where  c  is  a  function  of  the  determinant  of  the  spectral  density  matrix  Q,  which  is  defined 
to  be 


Q{kcoA)  =  TE 

'  G 
G 


Xx  (k)Xf(k)  Xi  (k)X$(k) 

_  x2(k)x;(k)  x2(k)x*(k) 

X\X\  {kuj  A)  G  X\X2  (k(jjA) 

GX2X2(ku!A ) 


defining  =  2n/T.  For  large  T,  Q(kuA)  — >  Q{f). 


The  inverse  is  easily  found,  with  the  determinant  given  by 

\Q(f)\  =  GX\x\  (. f)G  X2X2  (/)  -  \G  X±X2  (/)  I2 


=  G 


X\Xi 


(f)GX2X2(f) 


1  - 


\GXlX2(f)f 


GX \X\  (f)GX2X2(f) 


If  we  define  a  correlation  coefficient 


then  we  have 


712  (/) 


G 


X\X2 


(/) 


[GXlXl(f)G 

X2X2  C f )]* 


|Q(/)|=G;ci;ci(/)Ga;2;C2(/)(l-|7l2(/)|2) 


(B8) 

(B9) 


For  large  observation  time  T,  we  are  able  to  replace  the  Xi(k )  with  the  Fourier  transform 
Xi(2itk/T)/T  — >  Xif) /T  and  the  conditional  pdf  becomes 

p{X\Q)  =ce-12f**'(f)Q-1(f)x(f)  (B10) 


Taking  the  log,  and  considering  only  those  terms  dependent  on  the  parameter  td,  it  can  be 
found  that  the  maximum  likelihood  estimate  of  td  maximises  (using  the  notation  of  [19]) 

-j3 =2tJ  c„.j(/)|GU  (/)|  t  df  «»“) 
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where 

GxlX2(f)  =  ^X(f)X*(f)  (B12) 

is  the  estimate  of  the  cross  power  spectral  density  function  and  712  (/)  is  the  coherence 
[20,  21].  In  deriving  the  above,  we  make  use  of  the  relation 

GxlX2(f)  =  \GxlX2(f)\e-^ftd  (B13) 

It  is  straight  forward  to  derive  the  Cramer-Rao  bound  of  the  variance  of  td  [19]. 

The  log  likelihood  function  of  interest,  reduces  to 

m  =  - f  (B14) 

Taking  derivatives,  and  noting  that  E[GXlX2(f)\  =  GXlX2(f),  we  have  the  Cramer-Rao 
bound 

(m5) 
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